Dark matter in the solar system III: The distribution function of WIMPs at the Earth 

from gravitational capture. 
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In this last paper in a series of three on weakly interacting massive particle (WIMP) dark matter 
in the solar system, we focus on WIMPs bound to the system by gravitationally scattering off of 
planets. We present simulations of WIMP orbits in a toy solar system consisting of only the Sun and 
Jupiter. As previous work suggested, we find that the density of gravitationally captured WIMPs 
at the Earth is small and largely insensitive to the details of elastic scattering in the Sun. However, 
we find that the density of gravitationally captured WIMPs may be affected by external Galactic 
gravitational fields. If such fields are unimportant, the density of gravitationally captured WIMPs 
at the Earth should be similar to the maximum density of WIMPs captured in the solar system 
by elastic scattering in the Sun. Using standard assumptions about the halo WIMP distribution 
function, we find that the gravitationally captured WIMPs contribute negligibly to direct detection 
event rates. While these WIMPs do dominate the annihilation rate of WIMPs in the Earth, the 
resulting event rate in neutrino telescopes is too low to be observed in next-generation neutrino 
telescopes. 

PACS numbers: 95.35.+d,96.25.De,95.85.Ry,96.60.Vg 



I. INTRODUCTION 
A. Dark Matter Detection in the Solar System 

A number of lines of evidence point to the existence 
of a significant amount of dark matter in the universe, 
although its identity is a mystery [e.g., 1, 2, and refer- 
ences therein]. Perhaps the most popular candidate for 
dark matter is a species of WIMP, both because such 
particles appear naturally in extensions to the Standard 
Model of particle physics and because such particles have 
astrophysically desired properties [3-5] . 

There are experiments underway to identify WIMPs 
in collider experiments (e.g., the LHC [6-8]) or by their 
annihilation productions throughout the Galaxy [9-11], 
although it will be challenging to conclusively determine 
WIMP properties. In the case of collider experiments, 
WIMPs cannot be directly observed, and mapping any 
anomalous interactions to a particular theory will be dif- 
ficult. Astrophysical detection of WIMPs is complicated 
by both the uncertainty in the distribution of dark mat- 
ter in the densest parts of the Galaxy as well as poorly 
characterized foregrounds [e.g., 1, 12-15, and references 
therein] . 

There are several experiments to look for WIMPs in 
the solar system. Direct detection experiments look for 
the tiny (~ 10 — 100 keV) recoils of nuclei struck by 
WIMPs. Current experiments have ~ 10 kg of fiducial 
target mass. The XENONIO [16, 17] and GDMS [18] 
experiments currently have the best constraints on the 
spin-independent WIMP-proton a^^ and spin-dependent 
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WIMP-neutron cr^^ WIMP-neutron elastic scattering 
cross sections, at the levels of a^^ ^ 4 x 10"*'* cm^ and 
crf^ < 10-38 cm2 for WIMP mass - 100 GeV. Up- 
coming experiments should gain a factor of ^ 100 in sen- 
sitivity [19-22]. 

Neutrino telescopes are searching for neutrinos from 
WIMP annihilation in the cores of the Earth and the 
Sun. The current best constraint on the WIMP-proton 
elastic scattering cross section a^^ < 10^'^^ cm^ for 
m-^ ~ 100 GeV come from flux limits of neutrinos from 
the Sun [23, 24]. 

While the backgrounds in direct detection experiments 
and neutrino telescopes have been meticulously exam- 
ined, the astrophysical properties of WIMPs need to be 
understood either in order to make accurate predictions 
for the event rates in such experiments, or to constrain 
particle physics models from data. Thus, the distribu- 
tion function (DF) of WIMPs in the solar system needs 
to be characterized. Aside from uncertainties in the DF 
of halo WIMPs streaming through the [e.g. 25-27], there 
is theoretical uncertainty in the DF of WIMPs bound to 
the solar system. As has been demonstrated by several 
authors [28-31], this latter population may have a pro- 
found impact on predictions for direct detection event 
rates and on the annihilation rate of WIMPs captured in 
the Earth. 

In order to characterize the bound WIMP population, 
we have undertaken a program of simulating WIMP or- 
bits in the solar system, taking into account the possibil- 
ity of further scattering in the Sun. The results from our 
simulations in a toy solar system consisting of Jupiter on 
a circular orbit about the Sun are described in a series 
of three papers, of which this is the last. In Paper I [66], 
we simulated a population of WIMPs bound to the so- 
lar system by elastic scattering in the Sun, a population 
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originally postulated by Damour and Krauss [28]. We 
found that the DF of this population depended on both 
the WIMP mass and the strength of the WIMP-baryon 
interaction, but that the population was too small to 
significantly enhance direct detection event rates or to 
produce an observable signature of WIMP annihilation 
in the Earth. However, in Paper II of the series [67], 
wo found that the extended lifetime of WIMPs captured 
in the Sun had interesting consequences for searches for 
WIMP annihilation in the Sun. Both the gravitational 
interactions between WIMPs and planets and finite op- 
tical depth in the Sun to WIMPs altered the standard 
picture that all WIMPs captured in the Sun immedi- 
ately thermalize. These modifications to the standard 
thermalization pictured imply lower annihilation rates of 
WIMPs in the Sun for m-^ > 1 TeV or for low elastic 
scattering cross sections. 



B. Gravitationally Captured WIMPs 

In this paper, Paper III, we present simulations of the 
orbits of another class of bound WIMPs: those captured 

in the solar system by gravitational interactions with the 
planets. Previous work on this population has relied on 
treating all WIMP-planet encounters as local, and with 
only a crude treatment of WIMP-baryon encounters in 
the Sun. In order to make comparisons with our simula- 
tions, we briefly outline the previous work on the gravi- 
tationally bound WIMP population. 

The first to estimate the size of the gravitationally 
bound WIMP population was Gould [32, 33]. Gould ap- 
proximated all encounters between planets and WIMPs 
as local; gravitational scatters do not change the WIMP 
speed with respect to the planet, but do change the ori- 
entation of the orbit. This translates to a change in ve- 
locity in the heliocentric frame, and hence, to changes 
in the WIMP energy and angular momentum. Using a 
random walk approach, Gould found the average time for 
the angle between the direction of WIMP velocity with 
respect to the planet and the direction of motion of the 
planet to change by order unity as a function of WIMP 
speed, the "angular diffusion" timescale. He interpreted 
this as the timescale on which the planets could move 
WIMPs in or out of a particular region of phase space, 
sinc;c! whether a WIMP is bound or unbound to the solar 
system depends on the direction of the WIMP velocity in 
a planet-centric frame with respect to the planet motion. 
Goidd then estimated the DF of WIMPs at the Earth 
using the following detailed balance approximation. 

Assuming that Galactic halo WIMPs can be treated 
as having a Maxwellian DF near the solar system, and 
that the Galactic dark matter halo is non-rotating in 
an inertial Galactocentric frame, the DF of low speed 
WIMPs (ones that may be captured gravitationally) is 
nearly constant, f{v) « /. Since the escape speed from 
the solar system at the position of the planets is small, 
the phase space density of halo WIMPs at each planet 



is also « /. Gould argued that if the angular diffusion 
timescale were less than the age of the solar system, the 
phase space density of bound orbits should be the same as 
the phase space density of the unbound halo WIMPs for 
a given WIMP speed in an inertial frame moving with the 
planet. The flow of WIMPs filling the bound phase space 
is countered by the flow of WIMPs becoming unbound to 
the solar system. The angular diffusion timescale associ- 
ated with Jupiter is of order Myr for any part of phase 
space accessible to Jupiter. In Gould's picture, the phase 
space corresponding to bound Jupiter-crossing WIMPs 
should have the same density as the phase space asso- 
ciated with unbound orbits. Furthermore, Gould found 
that the timescale associated with the Earth and Venus 
for speeds with respect to the Earth of u < 30 km 
was less than the age of the solar system. Thus, for such 
speeds, the WIMP phase space density should be the 
same for any orientation of the velocity vector in an in- 
ertial frame moving with the Earth (geocentric). Some 
parts of phase space for it > 30 km s^^ is empty in this 
picture, but the majority of the accessible phase space at 
those speeds corresponds to unbound or Jupiter-crossing 
orbits. Hence, the speed distribution of WIMPs at the 
Earth in a frame moving with the Earth should be identi- 
cal to the speed distribution of halo WIMPs in free space 
(outside the potential well of the Sun) for it < 30 km s~^. 
Gould found that the free space approximation was rea- 
sonable for larger speeds, too. 

Gould neglected WIMP-nucleus scattering in the Sun. 
To determine the importance of this effect, Lundberg and 
Edsjo [31] also treated WIMP-planet encounters as local, 
and solved a gravitational diffusion equation for WIMP 
orbits in a solar system consisting of Jupiter, the Earth, 
and Venus. The Sun was either treated as a point mass 
or as infinitely optically thick to WIMPs. The timescale 
for hitting the Sun was estimated using a small set of 
individual WIMP orbit simulations. They found that the 
DF for WIMPs if the Sun were infinitely optically thick 
to WIMPs was substantially smaller than if the Sun were 
a point mass. 



C. This Work 

In light of previous work on gravitationally captured 
WIMPs, there are several reasons to perform suites of 
WIMP orbit simulations. First, the work of Lundberg 
and Edsjo [31] suggests that scattering in the Sun is an 
important loss mechanism for bound WIMPs. It would 
be useful to understand the degree of depletion as a func- 
tion of the WIMP optical depth in the Sun. Second, 
both Gould and Lundberg & Edsjo treat WIMP-planet 
interactions as local. However, these treatments neglect 
long-range encounters, short-period interactions, and is 
fundamentally insensitive to resonances (although these 
are incorporated for a set of Earth-crossing WIMPs in 
Lundberg & Edsjo ), which are known to be important in 
determining the dynamics of the population of WIMPs 
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bound to the solar system by elastic scattering in the 
Sun (Paper I) and of minor bodies in the solar system 
[34-37]. Finally, since WIMPs tend to be captured on 
initially very loosely bound orbits, they may be affected 
by external gravitational fields, which are known to be 
important in shaping the Oort cloud [38, 39]. 

In this work, we present simulations the gravitational 
capture and evolution of WIMPs in the solar system. As 
in Paper I, we simulate orbits in a toy model solar system 
consisting of Jupiter on a circular orbit about the Sun. 
The reasons for using a simplified system are twofold, 
(i) Since the integration algorithm is new, it is useful 
to check its performance in a simple system. The toy 
system we use has a constant of motion, the Jacobi con- 
stant; its constancy throughout the simulations was an 
indication of the accuracy of the integration scheme, (ii) 
Since Jupiter is by far the largest planet in the solar sys- 
tem, it should dominate the dynamics of WIMPs. These 
simulations provide a solid basis for understanding the 
dynamics of WIMPs in more complicated systems, which 
we hope to simulate in the future. 

We describe the simulations in Section II, and present 
the resulting DFs in Section III, which we construct us- 
ing a method outlined in Appendix A. In that section, 
we also discuss the DF in context of solar depletion and 
Galactic gravitational fields. We show the direct detec- 
tion and neutrino telescope event rates from the gravi- 
tationally bound WIMP population in Sections IV and 
V. In each of those sections, we compare the event rates 
from gravitationally captured WIMPs to the population 
of WIMPs bound to the solar system by elastic scatter- 
ing in the Sun, the subject of Paper I. In Section VI, we 
discuss the our results in context of other work. The key 
points of this work are summarized in Section VII. 



II. SIMULATIONS 

Simulations were performed using the algorithm de- 
scribed in Paper I, which we briefly summarize here. 

For most of the WIMP path, we integrated the orbits 
using a symplectic integrator optimized for systems in 
which one body dominates the gravitational potential of 
the system, and for which the gravitational force does not 
deviate significantly from an inverse square law [40, 41]. 
Symplectic integration is desired for long-term orbit inte- 
grations because errors are oscillatory instead of growing 
with time. This particular symplectic integrator is effi- 
cient for integrating the highly eccentric orbits character- 
istic of WIMPs because it allow for variable time steps 
(short at perihelion, long near aphelion) . To achieve this 
in a symplectic way, it treats time t and the WIMP en- 
ergy — po as conjugate variables; each step in time At is 
related to a step in the new fictitious time coordinate As 
by 



coordinates, respectively. For the choice 

GMq 



5(r,P) = - 



$(r,i)' 



(2) 



At = 5(r,p,t)As, 



(1) 



where r and p are the WIMP position and momentum 



where $(r, t) is the gravitational potential, Preto and 
Tremaine [41] show that the integrator exactly traces the 
solution to the two-body problem with only a phase error. 
Since typically <E> w — GMq/ |r — r©!, where Tq is the 
position of the Sun, At oc |r — TqI, so that for fixed As, 
shorter time steps are taken at perihelion than aphelion. 

We integrate the eight-dimensional equations of mo- 
tion using a second-order leapfrog mapping using a fixed 
fictitious time step As = h. We use only a second order 
integrator because we are interested in the behavior of 
an ensemble of orbits, and not the precise orbits of indi- 
vidual WIMPs. Even for this low order mapping, there 
is no numerical precession of orbits, and errors in the 
Hamiltonian are oscillatory in nature. 

Although it would be ideal use the symplectic inte- 
grator with fixed fictitious time step h throughout the 
integrations, it is too time consuming to be practical. 
This is because the integrator is not optimized to han- 
dle potentials that deviate significantly from <i>(r,t) = 
—GMq/It — TqI, which is the case in the interior of the 
Sun or when WIMPs experience close encounters with 
planets. In those cases, h would need to be set pro- 
hibitively small in order to resolve those potentials. 

Instead, we treat passages through the Sun and close 
encounters with planets using alternate methods, which 
allows for h to be set to a reasonably large value. 
While this breaks the Hamiltonian flow of the symplec- 
tic scheme, we have taken care to insure that our meth- 
ods for treating the special cases minimize errors in the 
Hamiltonian. 

For the passages through the Sun, we exploit the fact 
that tidal forces from the planets are much smaller near 
the Sun than they are elsewhere in the orbit. We treat 
passages through the Sun as a two-body problem, and are 
able to map the coordinates of the WIMP as it enters the 
Sun to its coordinates upon exit. In Paper I, we show that 
this method does not induce additional numerical errors. 

We define a sphere (or "bubble" ) around each planet 
in which we allow another break to the symplectic al- 
gorithm described above. In this bubble, we still inte- 
grate the WIMP orbits using the symplectic algorithm, 
but with a new fictitious time step h' tuned to achieve 
a minimum accuracy criterion for the Hamiltonian. To 
find this fictitious time step, we first integrate the WIMP 
orbit with the fiducial h for the region outside the planet 
bubble and outside the Sun. If a minimum accuracy re- 
quirement AE/E = \pQ + E{Y,t)\/E{r,t) is met, then the 
orbit is allowed to continue. If the minimum accuracy re- 
quirement is not met, the trajectory through the bubble 
is integrated with a slightly smaller h' . This procedure is 
iterated until the minimum accuracy requirement is met 
or the error plateaus. We tune the bubble size, the mini- 
mum accuracy requirement and h to minimize the overall 
integration time while maintaining small oscillatory er- 
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rors in the Hamiltoiiian tliroughout the integration. In 
later sections, we describe our specific choices for these 
variables. 



A. (Astro) Physical Assumptions 

The Solar System: We perform simulations of WIMP 
orbits in a toy solar system consisting of Jupiter {%) on a 
circular orbit of a% = 5.203 AU the Sun (modeled using 
[42]). This system admits a constant of motion, the Ja- 
cobi constant, which is a useful check on the accuracy 
of the integration algorithm. For simplicity, we treat 
Jupiter as having a constant mass density; this is not a 
realistic representation of Jupiter's structure, but only a 
tiny percentage of simulated WIMPs ever go through the 
planet. WIMP-baryon encounters in Jupiter are also ne- 
glected for similar reasons; in addition, the optical depth 
of Jupiter to WIMPs is negligible compared to that in 
the Sun. Since the kinetic energy and speeds of solar nu- 
clei are small compared to those of the WIMPs, we treat 
solar nuclei as being at rest with respect to the Sun. 

Dark Matter: The scattering probability of WIMPs in 
the Sun is completely determined by the solar model, the 
WIMP mass m^, and the cross sections and a^^ . 
Since we suspected that the bound WIMP DF would 
not strongly depend on scattering in the Sun, we chose 
only one point in the WIMP mass-cross section param- 
eter space to use for the simulations, = 500 AMU, 
= 0, and = 10"''^ cm^. This point lies below the 
best limits on WIMP parameter space from direct detec- 
tion experiments using standard assumptions about the 
halo WIMP DF [16, 18]. However, in order to extrapolate 
our DFs to other points in WIMP parameter space, we 
kept track of the integrated optical depth of each WIMP 
as a function of time. 



B. Starting Conditions 

In deciding how to arrange the initial conditions, it is 
useful to think about the flux of dark matter particles 
into a sphere of radius R centered on the Sun. The flux 
for an isotropic distribution function / is 

F{R,v) = 4Trv^f{vs{R,v))x^v cosed cose (3) 

= 'KV^f{vs{R,v))dvd{cos'^e), (4) 

where 7r/2 < 9 < w is the angle between the velocity v 
and the position vector R for incoming particles, and Vg 
is the speed of the particle relative to the Sun but far 
outside its gravitational sphere of influence. We have in- 
voked Liouville's theorem to find the WIMP phase space 
density at an arbitrary distance from the Sun. The total 
number of particles going inward through this spherical 
shell per unit time is 



It is useful to express this rate in terms of the specific 
energy E and specific angular momentum J instead of v 
and cos^ 6. Given that 

E = i^;2+ci>Q(i?) (6) 

J = Rvsine, (7) 

We find 

TV = tt/ (y2E^ dEdJ'^. (8) 

Therefore, the number of particles going through any 
shell is independent of the radius of the shell for a given 
energy and angular momentum; this is to be expected 
since there is no loss of particles between shells. 

If we were to sample all particles that fiow in towards 
the Sun, we would sample the energy according to to 
f{^/2E) and the angular momentum to be uniform in 
J^. However, by restricting the range of incoming par- 
ticles that are sampled to those that could be scattered 
onto boimd orbits, we can speed up the calculation. 

To find the range of E for which particles might pos- 
sibly be gravitationally scattered by Jupiter onto bound 
orbits, it is useful to think of gravitational capture in 
the following way. In the frame of the planet, the parti- 
cle speed does not change during the encounter, but its 
direction with respect to the direction of motion of the 
planet does. If the particle has a velocity v with respect 
to the Sun before encountering Jupiter, it will have an 
initial speed with respect to Jupiter of u = v — V'v, where 
v\ is the velocity of Jupiter with respect to the Sun. Af- 
ter encountering Jupiter, the particle will have a velocity 
u' with respect to Jupiter and v' = u' -|- 'V\ with respect 
to the Sun. For particles that were barely unbound to the 
solar system to begin with, it takes only a tiny deflection 
of the orbit to bind it to the solar system. However, for 
particles with increasingly higher energy with respect to 
the Sun, it takes an ever greater deflection by Jupiter to 
bind the particle. 

In order to find an upper limit to the energy from which 
particles may be captured, consider the most extreme en- 
counter possible. This is the case of a particle that has a 
tiny impact parameter with respect to Jupiter, and which 
has its initial velocity aligned with Jupiter's direction of 
motion. Therefore, the particle's velocity with respect to 
Jupiter is 

u = v — v%, (9) 

where v = )v|. The particle will be deflected through 
180°, so that 

u' = -{v-v%) (10) 
v' = 2v% - V. (11) 

The requirement that the particle is bound to the solar 
system after the scatter is equivalent to the statement 



N{R) = ATT^R^v^f{vs{R, t;))dt;d(cos2 6*). (5) 



v'\ < V2v%. 



(12) 
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Therefore, 



or 



and so 



2v\<v< {2 + V2)v%, 



(2 - V2)v'k. <v< 2v%, 



GMq 

a% 



(l + v^) 



(13) 

(14) 

(15) 
(16) 



This corresponds to a speed outside the gravitational 
sphere of influence of the Sun of 



= 2(1 + ^2) 
= 41 km/s. 



1/2 



v% 



(17) 
(18) 



No WIMP with a speed far from the Sun that exceeds 
41 km s~^ with respect to the Sun can be gravitationally 
captured. 

In addition to hmiting the range of E that we sample, 
we can also speed up the calculation by constraining the 
range of sampled. This constraint is equivalent to 
specifying a range of orbital perihelia to probe, given 
that the perihelion rp is defined by 



1 



(19) 



such that the angular momentum for a given energy E 
and perihelion is described by 



J^{E,rp)=2rp{Erj, + GMQ). 



(20) 



The goal is to make the range of rp (and hence, ,P) large 
enough to encompass all orbits that might become bound 
to the solar system while keeping the range small enough 
so as not to waste computing resources by following un- 
necessary orbits. 

Wc divide the gravitational scattering simulation into 
two parts, each defined by a different range of energy and 
perihelion: the "Regular run" and the "High Perihelion 
run." The Regular run samples particle orbits with: 



Q<E< 4/50 = -(44 km/s)2. 



Tp < 10 AU. (21) 



The maximum perihelion of 10 AU was chosen to be large 
enough — twice the semi-major axis of Jupiter — so that 

this run would contain the vast majority of particles that 
are gravitationally captured. If the Regular run misses 
any bound orbits due to the limit on r^, those orbits 
should be found in the High Perihelion run, defined by 



E < 4/50 



10 < Tp < 20 AU. 



(22) 



If we were to sample E and according to the distri- 
bution of particle energy and angular momentum squared 



flowing in towards the Sun, Eq. (8), the sampling prob- 
ability would be: 



G{E, J2) oc < 



G [J2(£;,r™"), J2(i;,r™''^)) 
0, J2 e [j2(S„„,r™"), J2(£;,r7")) (23) 



or J2 e [J2(£;,r™"^), J2(£;, 



min\ ^ 
min T ' p J ^ 



in the range Emin < E < Emax and J'^{E,, 
J2 < J'^iEmax, r™"^), where r^;*"^ and r^^"" are the max- 
imum and minimum perihelia allowed in each run. These 
ranges describe the maximum extent of E and for any 
given run. This sampling probability is highest in the 
high energy, high angular momentum part of the range 
considered. However, we want to sample proportionally 
more low energy orbits in both the Regular and High 
Perihelion runs, since these are most easily captured. We 
sample 



(24) 



G{E) = f{V2E) 



in the range Emin < E < Emax, and uniformly sample 
J2(£;,r™") < J2 < J2(£,r;""^). 

We treat the halo WIMPs has having a Maxwellian DF 
in Galactocentric coordinates in the solar neighborhood. 



/h(x,Vfc,i) = 



(27ra2) 



(25) 



where Vh is the WIMP speed in Galactocentric coordi- 
nates, far outside the sphere of influence of the Sim. 

~ Px/"^x is local WIMP number density, where 
the Px ~ 0.3 GeV cm~^ [25]. We set the one-dimensional 
WIMP velocity dispersion a = Vq/V^, where Vq « 
220 km s~^ is the speed of the Local Standard of Rest 
[43]. In heliocentric coordinates, the DF is 



(26) 



We use the angle-average of Eq. (26) to set the initial 

conditions. 

Once a sample particle's orbital parameters E and 
are selected, its initial position is determined by ran- 
domly orienting the position vector to a point on a spher- 
ical shell with fixed radius R relative to the Sun. The 
initial speed vector is chosen to be oriented inward, with 
the angle 9 relative to the position vector determined by 
J^. The speed v is fixed by R and since J = Rva'mO. 
The azimuth of the velocity vector relative to the position 
vector is also randomly chosen. Thus, the initial position 
and velocity of the particle are completely determined. 



C. Coordinate System Choice 

In general, we prefer to use heliocentric coordinates 
when integrating the equations of motion. However, the 
symplectic integrator breaks down at large heliocentric 
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TABLE I: The initial integration conditions for the gravita- 
tional capture simulation as a function of initial speed v and 
Kepler perihelion rp. The values of h are in units of Rq^ yr. 



Initial spoi'd r Jan s 


r^, < 5 Ai: 


r„ > .-) AU 


V <\Q 


2 X 10"^ 


3 X 10~^ 


10 < < 20 


5 X 10"'^ 


1 X 10-'^ 


u > 20 


1 X 10"^ 


2 X 10"® 









distances. This is because the gravitational potential in 
heliocentric coordinates has the form 



$(r,i) = - 



GMq ^ 



GM, 
r - Til 



+ 



,(27) 



where I denotes a planet. At large heliocentric distances, 
the last term in Eq. (27), the indirect term, becomes 
large compared to the gravitational potential of the Sun. 
Because we want to maintain g{v,t) « [r — FqI even at 
large distances from the Sun, we switch to barycentric 
coordinates far from the Sun. We use a crossover ra- 
dius Vc = 53 AU, which is point at which the value of 
the indirect term from Jupiter is approximately 10% the 
gravitational potential of the Sun. We find that using 
substantially smaller values of Tc induces large energy er- 
rors due to frequent breaks to the symplectic integration 
scheme. 



D. Setting h 



TABLE II: Choices for the fictitious time step ft as a function 

of semi-major axis for the gravitational capture simulations. 
The semi-major axis refers to bound particles unless otherwise 
indicated. 



a range [AU] 


h[R^' yr] 


< 0.75 


10-* 


0.75 < a < 1.1 


7 X 10"^ 


1.1 < a < 1.6 


6 X 10"^ 


1.6 < a < 3.5 


2 X 10"^ 


3.5 < o < 6.2 


1.5 X 10"^ 


6.2 < a < 13 


7 X 10"® 


13 < a < 22 


10"® 


22 < a < 30 


7 X lO"'' 


30 < a < 45 


6 X lO"'' 


45 < a < 120 


5 X 10"'' 


120 < a < 200 


4 X 10"^ 


200 < a < 500 


3 X lO"'' 


a > 500 or unbound 


2 X lO"'' 



E. In the Sun 

Once a WIMP is within 0.1 AU of the Sun, we check 
if its perihelion will lie within 2Rq of the center of the 
Sun, where Rq is the radius of the Sun. If not, the sym- 
plectic integration continues without interruption. If the 
WIMP does go through that region, though, we use the 
two-body map to evolve the WIMP through the region 
near the Sun. If the WIMP goes through the Sun, we em- 
ploy Monte Carlo techniques to determine if the WIMP 
scatters on a solar nucleus. These techniques are further 
described in Appendix C of Paper I. 



We initially set h according to Table I. This is suf- 
ficient for a "first pass" through the solar system. For 
long-term integrations, in order to both control errors 
near Jupiter and to speed up integration if particles set- 
tle onto tighter orbits, we reset h after the particles pass 
through the Jupiter bubble, h is actually reset at the first 
aphelion after passing through the bubble since we have 
empirically determined that this is the point in the orbit 
at which a change of h causes the minimum error. We 
then set h according to Table II. Since the semi-major 
axis of an orbit can change substantially throughout the 
integration, it is useful to occasionally change h to match 
o, either to speed up the integration or improve accu- 
racy. We allow h to change after each passage through 
the Jupiter bubble up to 10 Myr; however, to control 
for errors caused by breaking the symplectic nature of 
the integrator repeatedly, we only allow h to be reset if 
the energy changes by more than 20% through a Jupiter 
bubble passage after that time. Again, h is always reset 
at aphelion. 



F. Near Jupiter 

We set the accuracy criterion to \^E/E\ < 5 x 10~^ at 
the point at which the WIMP exits the Jupiter bubble. 
The bubble size was set to 1%^ =2.1 AU for particles with 
semi-major axes a < 100 AU, and Iq^ = 3.7 AU for orbits 
with either a > 100 AU or that were unbound. 



G. Stopping Conditions 

There were three circumstances in which orbit integra- 
tions were terminated: if the WIMP became unbound to 
the solar system, the WIMP rescattered onto an orbit of 
a < 0.3 (never to cross the Earth's path again), or the 
lifetime of the WIMP reached = 4.5 Gyr the lifetime 
of the solar system. 

In Table III, we show how many WIMPs are simulated 
in each of the Regular and High Perihelion runs. Sim- 
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TABLE III: Gravitational scattering simulations 



Name 




Regular 


4.8212 X 10'' 


High Perihelion 


3.994 X 10^ 



ulations were performed using computational resources 
at Princeton University. Each run required ^ 10^ CPU- 
hours on dual-core 3.2 GHz processors. 

III. DISTRIBUTION FUNCTIONS 

We constructed the DFs using the procedure outhned 
in Appendix A. For the Regular run simulation, the to- 
tal (bound -I- unbound) DP was derived from a total of 
369084 crossings within = 0.001 AU of the Earth's 
orbit, a result of integrating w 4.8 x 10^ particles with 
initial conditions distributed as in Eq. (24). Of those 
particles, 322441 particles were bound to the solar sys- 
tem for at least a short time, and 1224 of those bound 
orbits went through the Sun at least once. However, not 
a single particle was elastically scattered in the Sun. Of 
the 322441 particles that were bound to the solar sys- 
tem, only 5856 ever crossed the Earth's path (i? = 1 AU, 
\A — ^c), of which 772 also went through the Sun. There- 
fore, while only a small fraction of the bound orbits in 
this simulation contributed to the distribution function 
at the Earth, a large fraction of Sun-penetrating particles 
did. 

In the High Perihelion run (10 AU < Tp < 20 AU), 
there were only 9473 intersections with the Earth's orbit. 
Of the nearly 4 X 10^ particle orbits simulated for the High 
Perihelion Run, 64559 became temporarily captured in 
the solar system, 335 contributed to the bound DP, and 
54 went through the Sun. As in the Regular run, none of 
the particles going through the Sun were scattered onto 
smaller orbits. 

A total of 70943 WIMP crossings from the two runs 
were used to build up the bound WIMP DP. 

The WIMP DP from the simulations, including both 
bound and unbound WIMPs, is presented in Fig. 1. The 
DP is displayed in terms of the geocentric speed v (the 
speed of WIMPs relative to the Earth in an inertial frame 
moving with the Earth) and divided through by the halo 
WIMP number density rt^. It is normalized such that 
the number density of WIMPs near the Earth is given by 
/ dvv^f{v). We sum the DFs from each the Regular and 
High Perihelion runs and add the analytic DP of WIMPs 
(using Liouville's theorem) from the halo that were not 
in the energy and angular momentum windows used to 
set up the initial conditions for the simulations. In this 
figure, we have also plotted the free space distribution 
function (the WIMP DP outside the sphere of influence of 
the Sun) and the DP of imbound WIMPs for comparison. 



Distribution Function in tlie Frame of tlie Eartli 




10 20 30 40 50 60 70 



V [l<m/s] 

FIG. 1: The total distribution function from the gravitational 
capture simulations compared against several theoretical dis- 
tribution functions. 



In this figure, we have also plotted the DP predicted 
by Gould's detailed balance argument, discussed in Sec- 
tion IB. This DP includes both unbound WIMPs and 
Jupiter-crossing bound WIMPs, assuming the latter have 
the same phase space density as the unbound WIMPs. 
This is the most direct comparison to make with previous 
estimates of the DP; Lundberg and Edsjo [31] estimate 
the WIMP DP in a solar system with more planets. 

The DF from the simulations is fairly well fit by the 
detailed balance DF at low speeds, but the fit is poor 
at higher speeds. Moreover, we find that the DF of pro- 
grade bound WIMPs (those circulating in the same sense 
as Jupiter; the component of angular momentum perpen- 
dicular to the reference plane, J^, is positive) is larger 
than the DF of retrograde WIMPs (J^ < 0). Those DFs 
should be identical in the detailed balance approxima- 
tion. We believe these discrepancies are due to a violation 
of the key assumption in the detailed balance argument, 
that the phase space density of bound WIMPs depends 
on a single timescale, the angular diffusion timescale. 

Instead, we find evidence that the timescale for ejec- 
tion of WIMPs from the solar system is different from and 
shorter than the angular momentum diffusion timescale. 
If only one timescale governed energy and angular mo- 
mentum diffusion, we would expect that the distribution 
of the initial phase space coordinates of the WIMPs that 
built the bound DF at the Earth would be similar to 
the distribution of all bound WIMPs. As an example, 
we would have expected the distribution of the initial 
WIMP angular momenta for all WIMPs bound to the so- 
lar system to look similar to the distribution of the initial 
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FIG. 2: The initial angular momentum distributions for all 
bound WIMPs (upper two curves), and the distribution of the 
initial angular momenta for bound WIMPs that contribute to 
the WIMP DF at the Earth (lower two curves) . The distribu- 
tions are divided by whether the WIMPs were initially moving 
prograde or retrograde with respect to Jupiter. 



WIMP angular momenta for WIMPs that contribute to 
the DF at the Earth. In such a scenario, WIMPs with ini- 
tially high angular momentum would have enough time 
to lose enough angular momentum so that the perihelia 
of the WIMPs would lie within the Earth's orbit for a 
time before being kicked out of the solar system. 

We show these angular momentum distributions in 
Fig. 2, separated by whether the WIMPs were ini- 
tially prograde or retrograde. In the figure, we show 
,P = J'^/{2GMQa%), the square of the WIMP specific 
angular momentum divided by the square of the angular 
momentum for a WIMP traveling at the escape speed 
and reaching perihelion at Jupiter's orbit. We include 
all bound WIMPs in this plot, but we do not normal- 
ize the distribution to take into account the fact that we 
oversampled low energy orbits in the initial conditions 
(Section II B). If we had, the feature near = 1 would 
be more prominent, as it corresponds with WIMPs with 
small impact parameters with respect to Jupiter. Halo 
WIMPs with high initial energies must have small impact 
parameters on Jupiter in order to be captured to the so- 
lar system. However, the qualitative differences between 
the bound WIMP and Earth-crossing bound WIMP ini- 
tial angular momentum distributions are present in every 
energy interval. 

We find that distribution of the initial WIMP angular 
momenta of the WIMPs in the DF at the Earth is skewed 
towards small angular momenta relative to the distribu- 



tion of all bound WIMPs. The high angular momentum 
WIMPs cannot lose enough angular momentum to reach 
the Earth's orbit before they are ejected from the solar 
system. The effect is most pronounced for retrograde 
WIMPs. 

We considered that skew in the angular momentum 
distribution might be a result of the dependence of the 
WIMP lifetime in the solar system on the initial phase 
space coordinates. The reason for believing this might 
be a significant effect is that the cross section for WIMP- 
planet encounters is a function of WIMP speed with re- 
spect to the planet. Prograde WIMPs typically have 
small speeds with respect to the planet, and so the 
WIMP-planet cross section will typically be high. There- 
fore, the timescale to eject a WIMP will be short. For 
retrograde WIMPs, the relative speed of the WIMP in- 
creases as the angular momentum increases. The WIMP- 
planet cross section should be small, and the lifetimes 
should be longer. In general, WIMPs that are initially 
on prograde orbits will stay on prograde orbits, and like- 
wise for retrograde orbits; Jupiter simply cannot move a 
WIMP with large, positive onto an orbit with large, 
negative Jz- In the solar system as a whole, we find that 
the increase in the lifetime for retrograde WIMPs makes 
up for the smaller capture probability; there is an equal 
number of prograde and retrograde bound WIMPs in the 
solar system. 

However, the lifetime distribution of Earth-crossing 
WIMPs shows this not to be the case for the sample 
of Earth-crossing WIMPs. In Fig. 3, we plot the lifetime 
distributions of the WIMPs contributing the DF at the 
Earth as a function of the initial WIMP angular momen- 
tum. We show the lifetime distributions for four WIMP 
populations: those initially prograde with low angular 
momentum (defined as < 0.5 and Jz > 0); prograde 
with large angular momentum {,P > 0.5, Jz > 0); retro- 
grade with low angular momentum (J^ < 0.5, Jz < 0); 
and retrograde with high angular momentum {J^ > 0.5, 
Jz <0). We find that the median lifetimes for the WIMP 
populations are within a factor of three of each other, 
which is insufficient to explain the skewed angular mo- 
mentum distributions in Fig. 2. 

Using these figures, we can also explain why there is 
a difference between the DF of prograde and retrograde 
orbits at the Earth. In general, a WIMP that starts 
out on a prograde orbit stays prograde throughout its 
stay in the solar system, and a retrograde orbit stays 
retrograde. The prograde WIMPs are much more suc- 
cessful than retrograde WIMPs at reaching sufficiently 
low angular momenta such that their perihelia may lie 
inside the Earth's orbit for some time, as demonstrated 
in Fig. 2. The prograde WIMP DF is larger than the 
retrograde WIMP DF because the angular momentum 
diffusion timescale is less for the prograde WIMPs, even 
though the typical ejection time for a retrograde WIMP 
is only slightly longer than for a prograde WIMP. 

There is still the question of why the detailed balance 
DF fits the bound WIMP DF at small geocentric speeds. 
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FIG. 3: Lifetime distributions for initially prograde or retro- 
grade WIMPs. The solid lines show the lifetime distributions 
for WIMPs with P < 0.5, and dashed curves mark those with 
P > 0.5. 



To answer this question, we note that WIMPs populat- 
ing this part of phase space tend to have small semi- 
major axes and have angular momentum vectors nearly 
aligned with Jupiter's but small enough in magnitude to 
intersect the Earth's orbit. Using a crude diffusion argu- 
ment, one can show that the timescale for large changes 
to the WIMP angular momentum should be similar to 
ejection timescale for these WIMPs, but not for other 
orbits. Using an impulse approximation, the change to a 
WIMP's speed perpendicular to the direction of a planet 
in a planet-centric frame is 



Su 



GMf 
bu 



(28) 



where Mp is the planet mass, h is the impact param- 
eter, and u is the WIMP speed in the planet-centric 
frame. If the WIMP orbit is nearly radial in a helio- 
centric frame, which we generally expect for Jupiter- 
crossing WI MPs tha t have perihelia inside the Earth's 
T,, where v is the heliocentric WIMP 



orbit, u ~ 

speed and wp is the planet's circular speed about the Sun. 
Thus, the change to the heliocentric WIMP speed is of 
order 



5v 



GMp 

bv 



The change to the WIMP energy is of order 
6a 



5E_ 
E 



GMc 



-vdv, 



(29) 



(30) 



and the change to the angular momentum is 
5 J ^ apSv. 



(31) 



Using a random walk approximation, the change to 
either the energy or angular momentum (denoted as X 
below) goes as 



((AX)2) ^ 10N{SXy 



(32) 



where iV is the number of times a WIMP hits a planet 
in time i, and can be approximated by 



N . 



t 



(33) 



where is the WIMP orbital period. The factor of 10 in 
Eq. (32) includes the Coulomb logarithm, which we have 
otherwise ignored in this simplified random walk calcu- 
lation (for a more comprehensive treatment, see [44]). 
The rms change to the energy goes as 
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SO the ejection timescale goes as 



t,, - 0.1 



Mp 

Mo 



(34) 



(35) 



The rms change to the angular momentum goes as 

2 



{{5J)) ^ 10 



p ' 



(36) 



and the timescale for a WIMP of j ' 
pletely radial orbit is 



1 to reach a com- 



tj 



such that 



Px. 



ap 



(37) 



(38) 



Therefore, unless a ^ ap, the timescale for large changes 
in the angular momentum tj will be much longer than 
the ejection timescale. However, if a ^ ap, the angular 
momentum diffusion timescale will be approximately the 
same as the ejection timescale. Thus, the detailed bal- 
ance assumption that the energy and angular momentum 
diffusion timescales are the same and equivalent to the 
angular diffusion timescale is met, and the WIMP DF 
should resemble the detailed balance DF at the lowest 
geocentric speeds. 

In general, though, we find that while Jupiter is effi- 
cient at changing the energy of the WIMP orbits, it is 
quite inefficient at changing the WIMP perihelia. 

A similar phenomenon arises in another context in the 
solar system. It is thought that comets in the Oort Cloud 
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FIG. 4: (a) Distribution of times at which WIMPs cross the 
Earth's orbit. 



originate in the outer solar system, a = 4 — 40 AU. 
Through interactions with the outer planets, the energy 
of the objects is pumped up, such that a > 1000 AU. 
However, the perihelia stay nearly constant throughout 
this process. External gravitational fields (from pass- 
ing stars or molecular clouds) are required to move the 
perihelia of the comets outside the orbits of the planets 
[38, 39]. Again, we see the discrepancy between the ejec- 
tion timescale and the timescale to radically change the 
angular momentum of a body. 



A. Equilibrium Time 

In choosing to use the angle-averaged halo WIMP DF 
to set the initial conditions, we implicitly assumed that 
time equilibrium time for the bound WIMP DF (the time 
from the birth of the solar system beyond which the DF 
changes very little) was greater than the orbital period 
of the Sun about the Galactic Center, which is « 200 
Myr. The angle-averaged DF is approximately equal to 
the time-averaged DF due to the inclination of the plane 
of the solar system with respect to the Galactic plane 
(see Paper I and references therein). If the equilibrium 
time is shorter than the Sun's orbital period about the 
Galactic Center, the procedure to determine the DF, as 
described in Appendix A, would need to incorporate a 
treatment of the anisotropy in the halo WIMP DF (Eq. 
26). 

We show the distribution of times at which WIMPs 
cross the Earth's orbit {R — 1 AU, \z\ < Zc) in Fig. 
4. It is apparent that we have very few crossing times 
past ~ 200 Myr, which would imply that the equilibrium 



timescale for the bound WIMP DF should be less than 
the orbital period of the Sun about the Galactic Center. 
However, there is a possibility that the equilibrium time 
might be larger than it appears in Fig. 4. While there is 
a sharp peak in the crossing time distribution near ~ 1 
Myr, which is due to typical chaotic orbits, there some 
structure in the distribution of crossing times beyond 
10 Myr, and it is important to understand what types of 
orbits create this structure. 

Most of the Earth-orbit crossings beyond ~ 10 Myr 
are due to WIMPs that are temporarily stuck near mean- 
motion resonances, with a minority of the orbit-crossings 
coming from WIMPs initially captured onto large-a or- 
bits that scatter deep into the solar system at late times. 
In addition, these WIMPs show signs of also being on 
Kozai cycles. Kozai cycles are either librating or circulat- 
ing solutions about a type of secular resonance in which 
the rate of perihelion precession lo is small. The charac- 
teristic behavior of such orbits includes large swings in 
eccentricity and inclination of an orbit while the semi- 
major axis remains roughly fixed. The phenomenon of 
chaotic trajectories mimicking regular orbits near reso- 
nances for long times has been found in a number of sys- 
tems, most of them two-dimensional [36, 45, 46]. In the 
context of the solar system, such resonance-sticking has 
been found in simulations of Kuiper Belt objects (KBOs; 
[36]) and comets [46]. The peak in the High Perihelion 
crossing time distribution near 50 Myr is due to a single 
resonance-sticking particle. 

Given that there are only five resonance-sticking 
WIMPs that contribute significantly to the WIMP DF 
at the Earth past 10 Myr, and that only one contributes 
past 100 Myr, it is clear that this long-lifetime tail in the 
DF is poorly sampled. Since these few WIMPs account 
for ~ 15% of the bound WIMP Earth orbit crossings, 
it is important to understand how big (or small) the 
contribution of the resonance-sticking particles can be. 
There are two important issues in estimating the possi- 
ble size of the resonance-sticking DF: whether the orbits 
of the resonance-sticking WIMPs are typical of the pop- 
ulation as a whole, and what the lifetime distribution of 
the WIMPs is. The former can only be determined by 
more simulations. There are perhaps some insights from 
previous work into the latter point. 

The DF of resonance-sticking WIMPs beyond ~ 100 
Myr can be estimated in the following way. The rate at 
which a WIMP crosses the Earth's orbit can be described 
by Nc, which should be constant as long as the WIMP 
is stuck to a resonance. If the lifetime distribution is 
iV(> t) (X then the total DF beyond a time t can be 
estimated by 



fres{>t) cx I N{t')Nc{t')dt' (39) 

a>l 

cx <(log(toA), «=1, (40) 
a > 1. 
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Duncan and Levison [36], in their simulations of KBOs in 
a solar system consisting of the four outer planets, and 
with orbits of KBOs restricted to the plane, find that 
a = 1, in which case the resonance-sticking DF fres{> 
100 Myr) should be a factor of several greater than what 
was found in our simulations. Even in this case, the 
resonance-sticking WIMP DF will not be greater than 
the bound WIMP DF for typical chaotic WIMPs. Pre- 
vious work on planar systems, however, shows that usu- 
ally a > 1 [45, 46]. In that case, the only way we have 
underestimated the DF is if resonance-sticking WIMPs 
generically have a higher Nc than the WIMPs in our 
simulations. We note that previous work on resonance- 
sticking has focused on two-dimensional systems, so any 
extrapolation to fully three-dimensional systems should 
be treated with caution. 

For the gravitationally captured particles, the time av- 
eraging of the halo DF will turn out not to be justified 
if a > 1. If a < 1, the long-hfetime WIMPs will skew 
the equilibrium time higher, and so perhaps the time- 
averaging of the halo DF will be valid. In the former case, 
while the results and interpretation here are qualitatively 
correct, to make a precise prediction of the distribution 
function of gravitationally bound particles in the solar 
system, one should use the original, anisotropic halo dis- 
tribution function (Eq. 26) to translate the WIMP initial 
conditions in the simulations to a DF. 



B. Loss Mechanisms 

There are two means by which particles may be lost 
to the solar system other than gravitational scatter by 
planets: 

• Interactions with nuclei in the Sun (or, very rarely, 
the planets). Even though no bound orbits were 
scattered in the Sun in any of the gravitational cap- 
ture simulations, it is important to determine how 
the DF changes as a function of the strength of the 
dark matter-baryon interaction. 

• Interactions with external gravitational fields. 
Galactic tides and encounters with distant stars be- 
come important for bound orbits with a ^ 1000 
AU. Such Galactic gravitational fields are thought 
to be important in forming the Oort cloud as well as 
scattering Oort cloud comets into the solar system 
[38, 39]. It is important to understand how exter- 
nal fields will affect the distribution and lifetimes 
of WIMPs. 

In order to estimate the effect of the WIMP-nucleon 
cross section on the bound WIMP DF, we recorded the 
integrated optical depth as a function of time for each 
particle's orbit through the solar system. Very few of the 
bound orbits (1224/322441) ever went through the Sun, 
but the optical depths t of those that did are represented 
in Fig. 5. The median optical depth of Sun-crossing 
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FIG. 5: The distribution of total optical depth per particle 
of particles that enter the Sun. The solid line indicates the 
distribution for the Regular simulation, while the dashed line 
indicates that of the High Perihelion run. 



particles in the Regular run is Tmed ~ 10'^, and Tmed ~ 
2 X 10^"* in the High Perihelion run. Fig. 5 illustrates 
that very few particles have even a moderately high total 
optical depth if = 500 AMU, a^^ = 10"^^ cm^, and 
a^^ = 0. The WIMP-nucleon cross section (and hence, 
solar opacity) would need to be much, much higher in 
order for scattering in the Sun to rescatter any of the 
particles that pass through the Sun. 

To determine the maximum effect of scattering in the 
Sun, we found the bound DF of only those particles that 
never enter the Sun. This is represented by the trian- 
gles with error bars in Fig. 6. The majority of the 
bound WIMP DF is built up by particles that never en- 
ter the Sun. Therefore, the DF of bound WIMPs at the 
Earth depends only weakly on the strength or type of the 
WIMP-baryon interaction. 

The effects of the external gravitational fields are in- 
dependent of the WIMP mass and the WIMP-nucleon 
cross section. In order to estimate the consequences of 
these these forces, we assumed that Galactic tides pull 
the perihelia of all orbits crossing outward through 1000 
AU out of the solar system. This is approximately the 
radius at which the timescale for external fields to re- 
move orbital perihelia from the solar system is the same 
as the timescale for the planets to eject bodies [39]. In 
Fig. 6, we show the DF arising from particle-Earth orbit 
intersections that occur before the particle passes out- 
ward through 1000 AU (circles). The density of particles 
is noticeably lower that the total bound DF, generally 
by a factor of 3. It appears that WIMPs contributing 
to the DF at the Earth are initially captured on wide 
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FIG. 6: The DF of bound particles from the gravitational 
capture simulations. The squares mark the DF of all grav- 
itationally bound WIMPs. The circles indicate the DF of 
particles before they are lost to the solar system by Galactic 
tides. The triangles indicate the DF of particles that never 
enter the Sun. 



orbits that then shrink due to repeated encounters with 
Jupiter. Even though our treatment of external gravi- 
tational fields is crude, the DFs in Fig. 6 indicate that 
torques from the Galactic tide should be included in es- 
timates of the bound WIMP population at the Earth. 



C. Summary 

The main results of these simulations are twofold. 
First, the phase space density of unbound orbits is still 
quite a bit higher than that of the bound orbits above 
geocentric speeds w ^ 15 km s~^. We expect that this 
will be true even if anisotropic initial conditions are used, 
the external Galactic gravitational potential is more ac- 
curately modeled, and once better statistics of bound 
orbits are obtained. Second, the detailed balance DF 
is a poor fit to the WIMP DF for geocentric speeds 
w 30 km s~^ due primarily to the difference in ejec- 
tion and angular diffusion timescales, and secondarily to 
the presence of resonance-sticking orbits (points on the 
DF with larger-than-average error bars). Third, the DF 
is largely insensitive to rescattering in the Sun. Fourth, 
our crude treatment of Galactic gravitational fields sug- 
gests that these fields may be important in shaping the 
bound WIMP DF as well as the Oort cloud. 

Lastly, the phase space density of particles bound to 
the solar system by gravitationally scattering on Jupiter 
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FIG. 7: Comparison of the geocentric bound distribution 
functions for = 500 AMU and a^' = 10""^ cm^. The 
dot-dashed line indicates the phase space density of un- 
bound orbits. The squares show the results from the grav- 
itational capture simulations, the circles are the Large Mass 
DF (cTp^ = 0), and the solid magenta line indicates the esti- 
mated maximum DF resulting from spin-dependent scattering 
in the Sun {a^° = 10"^*^ cm^). 

is generally higher than that of particles bound by elas- 
tic scattering in the Sun ("solar captured WIMPs") for 
geocentric speeds d < 30 km s~^ and f > 50 km s^^. 
This is demonstrated in Fig. 7, in which we show the 
results of the gravitational capture simulation as well as 
DFs from the solar capture simulations of Paper I. If the 
spin-dependent cross section is high (Cp^ ^ 10~^° cm^), 
then the bound DF for the elastically scattered particles 
may be higher than the gravitationally captured parti- 
cles for 30 km s~^ < w < 50 km s~^, especially if the 
gravitationally captured WIMP population is depleted 
by external forces. However, the solar captured WIMP 
DF will be smaller or of approximately the same size as 
the gravitational capture DF if spin-independent scatter- 
ing dominates in the Sun, or if (7^^ < 10^"*° cm^. 

IV. DIRECT DETECTION 

Direct detection experiments look for nuclear recoil of 
rare WIMP-baryon interactions in the experimental tar- 
get mass. The WIMP-nucleus scattering rate per kg of 
detector mass per unit recoil energy Q can be expressed 
as [cf. 3] 



where daA/dQ is the differential interaction cross section 
between a WIMP and a nucleus of mass tua and atomic 
number A, and v is the velocity of the dark matter par- 
ticle with respect to the experiment. The lower limit to 
the integral in Eq. (41) is set to 

V„^^n = [niAQ /2fl\y/^ , (42) 

the minimum WIMP speed that can yield a nuclear recoil 

Q- 

We focus on direct detection rates for spin-independent 
interactions, but the results of this section can be ap- 
plied qualitatively to spin-dependent interactions as well. 
There is another class of direct detection experiment that 
is directionally sensitive [47-51]. In principle, the bound 
WIMPs should leave a unique signal in such experiments, 
but it would be challenging to distinguish this from the 
halo WIMPs given the small bound WIMP density, cur- 
rent errors in directional reconstruction, and high energy 
thresholds. 

We find direct detection rates assuming ^'^-'^Xe and ^^Ge 
targets, since the current and planned experiments most 
sensitive to the spin-independent (and spin-dependent 
neutron) cross section have multiple isotopes of either 
xenon or germanium as their target mass. We calculate 
the bound WIMP event rate for = 500 AMU and 
— lO^'''^ cm^. The event rate can simply be scaled 
for lower (or higher) spin-independent cross sections. The 
scaling for other values of and (7^^ is different, but 
can easily be determined. 

In Fig. 8, we show the differential direct detection 
event rate for the gravitationally bound WIMPs. For 
comparison, we also show the event rate predicted for 
halo WIMPs using an angle-averaged Maxwellian speed 
distribution. In addition, we show the direct detection 
rate for the detailed balance bound WIMP DF. As ex- 
pected, the event rate from the bound WIMPs in the 
simulation is somewhat less than predicted from detailed 
balance estimates. The maximum contribution of the 
bound WIMPs to the event rate is at Q = 0, at which 
point it is approximately ~ 0.3% that of the halo. In or- 
der to estimate the contribution of gravitationally bound 
WIMPs to the event rate in current experiments, we 
show the analysis windows for the XENONIO (shaded 
region) and CDMS (right of the vertical dashed line) ex- 
periments, which have xenon and germanium targets, re- 
spectively. The CDMS experiment should be completely 
unaffected by bound WIMPs; Vmin is larger than the 
maximum bound WIMP speed at the analysis threshold. 
The XENONIO experiment should be sensitive to bound 
WIMPs; however, the contribution to the total event rate 
is negligible. 

In Fig. 9, we compare the direct detection rate of 
gravitationally bound WIMPs to the total bound WIMP 
event rate, which includes the maximum contribution to 
the event rate from the solar captured WIMPs discussed 
in Paper I. The largest solar captured WIMP DF occurs 
for ~ a few hundred GeV, and if the spin-dependent 
WIMP-proton cross section afP > 10"^" cm^ and dom- 
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FIG. 8: The differential direct detection signal from the halo, 
the gravitationally bound WIMP population, and the detailed 
balance estimate for the toy model solar system. The shaded 
region indicates the XENONIO analysis region [16], and the 
vertical dashed line indicates the lower limit to the CDMS 
analysis window (which extends to Q = 100 keV) [18]. 
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FIG. 9: The maximum contribution to the differential direct 
detection rate for = 500 AMU and = 10""^ cm^. The 
upper lines represent the spin-independent event rate of halo 
WIMPs assuming a ^^^Xe and '^^Ge) targets. The lower solid 
lines show the event rate of gravitationally bound WIMPs; 
the lower dashed lines also include WIMPs bound to the solar 
system by solar capture assuming (t^^ = 10"'^'' cm'^). 
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inates the optical depth to WIMPs in the Sun. In that 
case, we find that the bound WIMP event rate is dom- 
inated by solar captured WIMPs; the maximum value 
of the differential event rate also occurs at Q = 0, and 
is ~ 0.5% the halo event rate. Since the solar captured 
DF is large at 30 < u < 50 km s~^ and small for other 
geocentric speeds, the gravitationally captured WIMPs 
dominate for Q 5 keV, which means they will domi- 
nate the bound WIMP signal in the XENONIO analysis 
window. 

However, the solar captured WIMP DF is typically 
smaller unless > 10"*" cm^ or > 10~^^ cm^. 

For smaller WIMP-proton cross sections, gravitationally 
captured WIMPs will dominate the bound WIMP direct 
detection signal unless the Galactic tidal fields or the 
gravitational potentials are strong enough to severely re- 
duce the gravitationally bound WIMP DF, as shown in 
Fig. 6. 

The main conclusion in this section is that the con- 
tribution of bound WIMPs to the event rate expected 
in direct detection experiments is negligible and will not 
affect parameter estimation based on the shape or nor- 
malization of the direct detection event rate. 



speed after scattering Vf must be less than the local es- 
cape velocity Vesc{^)- The second term in Eq. (43) is 
twice the annihilation rate 



r = (av) 



(45) 
(46) 



where {av)a is the velocity-averaged annihilation cross 
section, n{r,t) = Nn{r,t) (J cfrh = 1) is the density 
of WIMPs in the Earth. The factor of two in Eq. (43) 
comes from the fact that most popular WIMP candidates 
are self-annihilating. The coefficient is a constant as 
long as n(r, t) is time-independent. 

If the capture rate C is time independent (i.e., the dis- 
tribution function is time- independent), the annihilation 
rate can be calculated analytically: 



-Ctanh^(t/te), 



where 



te = (CCa) 



-1/2 



(47) 



(48) 



V. NEUTRINOS FROM WIMP ANNIHILATION 
IN THE EARTH 

WIMPs may accumulate and annihilate in the Earth. 
The signature of WIMP annihilation will be GeV to TeV 
muon neutrinos, which may be observed in terrestrial 
neutrino observatories (e.g., Antares [52], IceCube [53]) 
via the Cerenkov radiation of muons created in charged- 
current interactions of muon neutrinos in and around the 
experiment. In this section, we estimate the range of 
possible event rates due to WIMPs bound to the solar 
system assuming a neutralino WIMP. 

The muon flux in the telescopes is proportional to the 
annihilation rate F of WIMPs in the Earth. If WIMPs 
quickly settle into an equilibrium distribution in the 
Earth once they are captured, the annihilation rate may 
be found by solving 



N = C-2T, 



(43) 



where N is the number of WIMPs in the Earth. The 
capture rate of WIMPs in the Earth by elastic scattering 
is defined as 

C = /d^x / d^vdn V ^n^(x)^ 

X /(x,v,i). (44) 

Here, dcr^/dO is the WIMF-nuclcus clastic scattering 
cross section for nuclear species A and v is the relative 
speed between the WIMP and a nucleus. The number 
density of species A is described by n^(x). The cutoff 
in the velocity integral reflects the fact that the WIMP's 



is the equilibrium timescalc. In the limit that the equi- 
librium timescale is small or large relative to the age of 
the solar system t©, 



fft0Ae»l 
XlC'Catl fft0Ae<l- 



(49) 



For parts of WIMP phase space not experimentally ex- 
cluded, tQ, so the annihilation rate of WIMPs in the 
Earth is a sensitive function of the capture rate. Since 
the Earth has a shallow potential well (the escape speed 
from the center of the Earth is Vesc ~ 15 km s~^), only 
low speed WIMPs may be captured by the Earth unless 
the WIMP mass is close to the mass of one of the dom- 
inant nuclear species in the Earth [32]. Therefore, even 
though the bound WIMP population is small, it domi- 
nates the event rate in neutrino telescopes. This is illus- 
trated in Fig. 10, in which we show capture rates assum- 
ing CTp = 10~^^ cm^ using the Earth models in Refs. En- 
cyclopaedia Britannica [54], McDonough [55]. The cap- 
ture rate from only unbound halo WIMPs in shown with 
the solid line, which drops to zero for m.^ > 400 GeV 
due to the cutoff in the DF at the escape speed from the 
solar system. For reference, wc have also plotted the cap- 
ture rate for the free space DF, which the DF preferred by 
Gould [33] assuming gravitational diffusion from all plan- 
ets in the solar system. This capture rate is substantially 
larger for > 70 GeV than that found in our simu- 
lations because of the relatively large free space WIMP 
phase space density at small geocentric speeds. The max- 
imum capture rate from solar captured WIMPs (with 
(Tp^ consistent with supersymmetry values) in addition 
to the unbound WIMPs in shown with the long-dashed 
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FIG. 10: Capture rate of WIMPs in the Earth as a function 
of WIMP mass for a^' = 10"*^ cm^ 



line. This capture rate is smaller than that for grav- 
itationally captured WIMPs (dot-dashed Hne) because 
the gravitationally captured WIMP DF extends to lower 
speeds than the solar captured WIMP DF, and is gener- 
ically larger for geocentric WIMP speeds v < 30 km s^^ 
(see Fig. 7). The gravitationally captured WIMPs dom- 
inate the overall capture rate in the Earth for ^100 
GeV. We do not show the capture rate predicted by de- 
tailed balance because the DF is virtually indistinguish- 
able from the simulation WIMP DF at the small speeds 
relevant for capture in the Earth. 

To estimate a plausible range of muon event rates given 
the capture rates in Fig. 10, we explore a hypersurface 
of minimal supersymmetric model (MSSM) parameter 
space assuming the WIMP is a neutralino. We can in 
principle explore other models, but the MSSM yields, 
on average, somewhat larger spin-independent cross sec- 
tions. Given that iron is the most common element in 
the core of the Earth, and oxygen, silicon, and magne- 
sium the most common element in the mantle, none of 
which has spin-dependent interactions with WIMPs, only 
in WIMP models with appreciable spin-independent in- 
teractions will capture in the Earth be relevant. 

To estimate the neutrino-induced muon event rate for 
neutrino telescopes from neutralino annihilation in the 
Earth, we use routines from the publicly available Dark- 
SUSY v.5.0.2 code [56]. Because searching the space of 
the large number of free parameters in the MSSM is a 
nearly impossible task, DarkSUSY has a simplified set 
of inputs from which all other MSSM parameters are set 
in a physically motivated way. The seven free param- 
eters, specified at the weak-breaking scale, are: ^, the 



Higgsino mass parameter; tan/3, the ratio of the Higgs 
vacuum expectation values; M2, the mass of one of the 
gauginos, through which the other two gaugino masses 
are specified; mcp, the mass of the CP-odd Higgs (usu- 
ally denoted by rriji, which we avoid in order to prevent 
confusion with m^, the mass of a nucleus with atomic 
number A); toq: which sets the masses of the lepton and 
quark superpartners; and At and Af,, which parametrize 
the strengths of the trilinear couplings in the most gen- 
eral MSSM Lagrangian. 

To generate a set of MSSM models for the neutralino, 
we scan a seven-dimensional hypersurface of the MSSM. 
The range used for each parameter is given in Table IV. 
For jj,, M2, mcp, tan/3, and toq, we sample the range 
logarithmically, and sample the other parameters linearly 
in their ranges. We accept a model if it makes it through 
the colhder constraints, 0.05 < Vl^h^ < 0.125, and cr^^ > 
10^^^ cm'^. The upper limit on the allowed region of 
il^ft,^ is approximately the 3cr range of Q.dmh'^ from the 
WMAP-b analysis [2]. The lower limit is about half the 
3cr lower limit from that analysis, since the neutralino 
may not be the only dark matter species. We use 780 
models which satisfied the requirements in our scans for 
the discussion below. 

To estimate the muon event rate in a neutrino tele- 
scope, we set the muon energy threshold to E*-^ — 1 
GeV. This is somewhat optimistic for the IceCube exper- 
iment [31, 57] unless muon trajectories lie near and ex- 
actly parallel to the PMT strings, but it is reasonable for 
the more densely packed water experiments (e.g., Super- 
Kamiokande). We assume that the material surround- 
ing the detector volume, the target material for neutrino 
interactions, is either water or ice, since the largest cur- 
rent and upcoming neutrino telescopes are immersed in 
oceans or the Antarctic ice cap. We include all muons 
oriented within a 30° cone relative to the direction of the 
center of the Earth. 

We present muon event rates in neutrino telescopes 
for various DFs in Fig. 11. We show the event rates 
for WIMPs unbound to the solar system, as well as 
for gravitationally captured WIMPs (in addition to the 
halo WIMPs) and for the maximum DF from solar cap- 
ture (in addition to the halo and gravitationally cap- 
tured WIMPs). In the last case, we use the DF when 
(Tp'° = 1.3 X 10"'^^ cm^, which is nearly the maxi- 
mum spin-dependent cross section found in our param- 
eter scans. The solid black line in this figure represents 
the most optimistic flux threshold for IceCube [31, and 
references therein] . To show how the event rates depend 
on the SUSY models for a given spin-independent cross 
section, we mark the models on the figure according to 
which direct detection experiments bracket the cross sec- 
tion for a given neutralino mass. The open circles corre- 
spond to SUSY models with a^^ above the that lie above 
the 2006 CDMS hmit [58] . The triangles are models for 
which a^' hes between the 2006 CDMS limit and the 
current best limits on (a combination of XENONIO 
[16] and CDMS [18] limits), and squares denote models 
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TABLE IV: Ranges of parameters for the seven-parameter DarkSUSY MSSM inputs at the weak scale 



SUSY parameters: 


11 [GeV] 


Ma [GcV] 


mcp [GeV] 


tan/3 


mo [GeV] 


At 


At, 


min 


-50000 


-50000 


1 


1 


50 


-3 


-3 


mcix 


50000 


50000 


50000 


60 


20000 


3 


3 



consistent with all current direct detection experiments. 

While we find that the bound WIMPs, especially those 
gravitationally captured to the solar system, do increase 
in the muon event rate in neutrino telescopes, the event 
rates fall far below threshold for the models in our scans. 
While we cannot say that it is impossible for neutralino 
WIMPs to be observed by IceCube or other km^-scale 
experiments, since we are only sampling a small part of 
the SUSY parameter space, the prospects do not look 
good. If Galactic gravitational fields are important in 
the solar system for heliocentric distances as small as 
^ 1000 AU, the event rate due to bound WIMPs would 
be a factor of ~ 10 smaller yet, since the DF of WIMPs 
before crossing outward through r = 1000 AU is a factor 
of ^ 3 smaller than in the absence of such fields for low 
geocentric speeds, and P oc for such small capture 
rates. 

As a caveat, only the flux of muons created outside 
the detector volume is calculated in DarkSUSY. This 
was historically done because muon path lengths were 
long compared to the detector dimensions. More recently. 
Bergstrom et al. [59] found that muons created inside the 
detector volume dominate the signal for smaller WIMP 
masses (to^ ^ 300 GeV) in large (km'^) telescopes, and 
that the expected event rate from muons created within 
the detector volume depends quite sensitively on the con- 
figuration of detectors inside that volume. Therefore, 
the event rates used here ought to be considered a lower 
limit to the actual event rate in a large detector for neu- 
tralino masses < 300 GeV. From Fig. 11, we find 
that the event rate of muons created in the telescope vol- 
ume would need to be at least two orders of magnitude 
larger than the event rate of muons created outside the 
telescope in order to be observed. The ratio between 
the gravitationally bound WIMP event rate and the halo 
event rate, however, would be unchanged in this scenario. 



VI. DISCUSSION 

A. Comparison with Previous Work 

In this section, we discuss oxu results in comparison 
with previous work on gravitationally captured WIMPs, 
namely the work of Gould [33] and Lundberg and Edsjo 
[31]. First, we focus on direct comparisons between our 
simulations and the work on the toy solar system. Sec- 
ond, since the other authors considered the effects of 



planets, which we have not yet discussed in the context 
of the simulations, we describe what results from the toy 
solar system wc expect will hold in the true solar system, 
and what we expect might change. 

Toy solar system: First, we describe which of our re- 
sults agree with previous work, and then how they dis- 
agree. Our results agree with previous work in two im- 
portant ways, (i) While Lundberg & Edsjo find that 
scattering in the Sun may drastically reduce the WIMP 
DF in a three-planet (Jupiter, Earth, Venus) solar sys- 
tem relative to the DF if the Sun were a point mass, they 
find that the WIMP DF in a one-planet (Jupiter) solar 
system is largely unaffected by the details of scattering in 
the Sun. Our results (Fig. 6) confirm their findings, (ii) 
Even though our bound WIMP DF is smaller than the 
detailed balance DF, the two DFs are nearly identical for 
small geocentric WIMP speeds. Since the capture rate 
of WIMPs in the Earth is most sensitive to the density 
of the lowest speed WIMPs, our predictions for the event 
rate for neutrino telescopes match those predicted using 
detailed balance arguments. 

We find several important deviations from the detailed 
balance picture of Gould. First, we showed that angu- 
lar momentum and energy diffuse at different rates in 
the solar system. This implied a deficit in the bound 
WIMP DF for large geocentric WIMP speeds, as well 
as an asymmetry between the prograde and retrograde 
WIMP DFs. 

Second, we found a set of long-lived {t > 10 Myr) 
resonance-sticking WIMPs that contributed ~ 15% of the 
bound WIMP DF. Since our sample was small, we were 
unable to determine how statistically significant that con- 
tribution was. The resonance-sticking bound WIMP DF 
depends on (i) the orbital properties of the WIMP while 
stuck to a resonance, (ii) the distribution of WIMPs 
among resonances, and (iii) the lifetime distribution. The 
former two points will likely only be addressed in future, 
larger simulations. We used the lifetime distributions 
from studies of comets in the solar system to argue that 
the resonance-sticking WIMP DF would be at most a 
similar size to the total gravitationally captured bound 
WIMP DF found in this work. However, most of the work 
on comets was done in nearly planar systems, and with 
initial conditions for the comets that are quite different 
that the WIMP initial conditions [36, 46]. We caution 
that WIMP orbits in the solar system are fully three- 
dimensional, and that the WIMP DF depends not on the 
overall lifetime distribution of WIMPs in the solar sys- 
tem, but on the lifetime distribution of WIMPs on Earth- 
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FIG. 11: Muon event rates from (a) halo WIMPs unbound to the solar system, (b) halo and gravitationally bound WIMPs, 
and (c) halo and all bound WIMPs. Open circles mark MSSM models for which is above the 2006 CDMS limit [58] , filled 
triangles mark those with limits between that limit and the current best limits on cTp^ (set by XENONIO for < 40 GeV 
[16] and CDMS for > 40 GeV [18]), and filled squares denote models consistent with the best limits on elastic scattering 
cross sections. The solid line is an optimistic detection threshold for the IceCube experiment [31, and references therein]. 



crossing orbits. As an aside, all of the resonance-sticking 
WIMPs were originally captured onto orbits a > 500 
AU, such that none of these WIMPs would have con- 
tributed to the WIMP DP if Galactic tidal fields were 
strong at r > 1000 AU. The only way to determine how 
important resonance-sticking is in the solar system is to 
perform larger simulations than those presented in this 
work, and incorporating a better treatment of external 
gravitational fields. 

Finally, our work is the first to explore the possibility of 
Galactic gravitational fields affecting the bound WIMP 
DP. Even though our treatment of the fields is crude 
(removing WIMPs from the inner solar system as soon 
as they pass outward through r = 1000 AU from the Sun), 
it suggests that external gravitational fields may play a 
significant role in shaping the bound WIMP distribution 
in the solar system. Future simulations should include 
a more realistic treatment of the Galactic gravitational 
fields in order to make precision predictions for the bound 
WIMP DP at the Earth. 

The potential importance of the Galactic gravitational 
fields brings up a shortcoming of treating all WIMP- 
planet encounters as local. We find that the WIMP DF 
is much smaller when a crude treatment of Galactic tides 
is used because long-range capture of barely unbound 
WIMPs is quite important, and these WIMPs are typi- 
cally captured onto barely bound orbits. Therefore, the 
WIMP DF is sensitive to the details of capture in the 
solar system, which are not well-described by treating all 
encounters between WIMPs and planets as strictly local. 

The effect of more planets on the DF: Our conclusions 
are based on simulations in a toy solar system, while the 
true solar system is far more complex. The question is, 
how will putting Jupiter on a more realistic eccentric or- 



bit and the presence of other planets alter the DF? There 
are several things we expect. First, we expect the asym- 
metry between prograde and retrograde orbits to persist 
in a more realistic solar system since all planets revolve 
in the same sense about the Sun as Jupiter. Second, 
we expect the DF of resonance-sticking, Jupiter-crossing 
WIMPs of the types found in the toy solar system to 
be reduced, although it is not clear by how much. The 
DF of such WIMPs is expected to be smaller because en- 
counters with other planets (especially the outer planets) 
can perturb the WIMPs off the quasi-regular orbits they 
have when stuck to resonances. Although simulations of 
KBOs have shown that up to 1% of orbits can survive 
~ Gyr with a lifetime distribution A'^(> i) oc in a 
solar system consisting of only the four outer planets, 
we do not expect those results to generalize to WIMPs. 
Even in the toy solar system, only ^ 0.1% of WIMPs 
contributing to the DF were on resonance-sticking orbits 
lasting > 10 Myr. In addition, the typical WIMP peri- 
helia are much smaller than the KBO perihelia, so that 
WIMPs may experience close encounters with gas giants 
other than Neptune (the planet which governs much of 
the KBO dynamics). 

Finally, we expect the overall gravitationally bound 
WIMP number density to be within a factor of a few 
of what we found in the toy solar system, although this 
depends on whether a significant long lifetime tail (a re- 
sult of various resonances involving WIMPs interior to 
Jupiter's orbit) exists. The implications for direct de- 
tection experiments and neutrino telescopes will depend 
on how those WIMPs are distributed in velocity space. 
Here, we consider a few possible ways in which a more 
complicated solar system might affect WIMP orbits and 
the size of the WIMP population. 
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First, we consider the capture of WIMPs to the solar 
system. The gravitational cross section cx M|>, where 
Mp is the mass of planet P. The inner planets are not 
nearly massive enough to significantly boost the capture 
rate. In the outer solar system, Jupiter is the most mas- 
sive planet by a factor of ^ 3.3, implying that its gravita- 
tional cross section is at least a factor of ^ 10 higher than 
any other planet. Even if an outer planet were primarily 
responsible for capturing a WIMP to the solar system, 
there are two reasons why the WIMP DF at the Earth 
will be largely unaffected. First, if an outer planet cap- 
tures a WIMP on a Jupiter-crossing orbit, Jupiter will 
dominate the dynamics of the WIMP. Second, Jupiter is 
the innermost outer planet. Thus, the typical angular 
momentum of a WIMP captured by another planet will 
be higher than a WIMP captured by Jupiter. We de- 
termined in Section III that the diffusion of angular mo- 
mentum is typically much slower than energy diffusion. 
Even if a WIMP with initially large angular momentum 
diffuses down to Jupiter's orbit, Jupiter will again dom- 
inate the dynamics of the WIMP, and will tend to eject 
the WIMP before it can diffuse further in angular mo- 
mentum. Therefore, even if an outer planet captures 
WIMPs to the solar system, it is unlikely to affect the 
DF of WIMPs at the Earth. 

Next, we consider the combined effects of gravitational 
diffusion and rescattering in the Sun. While Gould ar- 
gued that the WIMP DF should be almost identical to 
the free space DF in the geocentric frame due to gravita- 
tional diffusion, Lundberg & Edsjo find that the DF is not 
much larger than the DF we found in our simulations. In 
order to incorporate the effects of scattering in the Sun in 
the WIMP diffusion equation, Lundberg & Edsjo started 
2000 WIMP orbits on a grid in geocentric velocity space, 
integrating orbits up to 49 Myr. For each WIMP that hit 
in the Sun within that time, the original point in veloc- 
ity space was assigned a scattering frequency {ly — 1/i/, 
where ti is the lifetime of the WIMP before scattering). 
The loss term in the diffusion equation is —nu, where n 
is the WIMP orbit density. This means that WIMPs are 
lost from the solar system on the timescale for them to 
hit the Sun. The treatment of scattering in Lundberg & 
Edsjo encompassed the various ways by which WIMPs 
are driven into the Sun. Simulations of near-Earth ob- 
ject (NEO) and asteroid orbits suggest that secular res- 
onances (occuring when either the rate of change of the 
longitude of perihelion vj or the longitude of the ascend- 
ing node Cl are the same as for one of the planets), as 
well as mean- motion and Kozai resonances, drive them 
into the Sun on ^ 1 — 10 Myr timescales [35, 37]. 

While these effects are incorporated into the diffusion 

equation of Lundberg & Edsjo, we find that there are sev- 
eral ways in which a full orbit integration with a Monte 
Carlo treatment of scattering could improve on the work 
of Lundberg & Edsjo, and reasons why the DF found by 
Lundberg & Edsjo is likely to be too small. First, as we 
found in Paper I and in Section III of this work, WIMPs 
can survive many passages through the Sun before be- 



ing removed from Earth-crossing orbits. Therefore, the 
lifetimes of WIMPs could be substantially longer that 
assumed by Lundberg & Edsjo. The exact amount by 
which the lifetimes would be extended depends on how 
many passages WIMPs make through the Sun each time 
they are driven into the Sun by the planets, and how 
deeply into the Sun the WIMPs penetrate. In addition, 
WIMP orbits may prcccss rapidly in the Sun, which may 
affect the type of orbit they are on. This may change the 
frequency with which WIMPs are driven into the Sun in 
the future. 

One effect that may decrease the DF is related to the 
short length of the Earth-crossing orbit simulations. For 
a large swath in the geocentric velocity space, WIMPs 
were neither ejected from the solar system nor driven 
into the Sun on timescales less than 49 Myr. However, 
WIMPs may be driven into the Sun on longer timescales. 
In particular, WIMPs on highly eccentric or highly in- 
clined orbits can survive a long time before being scat- 
tered onto orbits that lead either to ejection or penetra- 
tion of the Sun. This is because the orbits are almost per- 
pendicular to the direction of motion of the planets, lead- 
ing to high speed encounters with the planets in which 
the WIMPs are not strongly deflected. Therefore, Lund- 
berg & Edsjo may be underestimating scattering in the 
Sun by not simulating such orbits long enough. How- 
ever, the size of the effect will be determined by the op- 
tical depth in the Sun to WIMPs, and what the typical 
integrated optical depth is each time a high eccentric- 
ity or high inclination WIMP is driven into the Sun. If 
the integrated optical depth is small, then Lundberg & 
Edsjo may not be underestimating the effects of scat- 
tering for this population. Underestimating the lifetime 
of the WIMPs that encountered the Sun within 49 Myr 
will be a bigger effect than underestimating the scatter- 
ing probability of the WIMPs that neither are ejected 
nor pushed into the Sun on 49 Myr timescales since the 
latter is closer to the lifetime of the solar system. 

Finally, the distribution of times at which WIMPs en- 
ter the Sun for the first time is poorly sampled in Lund- 
berg & Edsjo , and the system is chaotic, so interpolat- 
ing the lifetimes among grid points in velocity space may 
not be the best way to interpolate lifetimes. In Paper 
I, we found that the DF of solar captured WIMPs at 
the Earth was dominated by the long lifetime tail in the 
Earth-crossing WIMP distribution. In fully integrating 
the WIMP orbits, the impact of the long lifetime, gravi- 
tationally captured WIMPs on the DF at the Earth can 
be understood. 

In summary, we suspect that the DF predicted by 
Lundberg & Edsjo is too conservative; for most of WIMP 
parameter space consistent with experimental limits on 
the WIMP-baryon cross section, the optical depth per 
passage through the Sun is small. Thus, WIMPs can 
survive many passages through the Sun before being 
removed from Earth-crossing orbits. The effect of the 
longer lifetimes to scattering on the DF will also de- 
pend on the ejection timescale. If the ejection timescale 
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is shorter than the scattering timescale, the WIMP DF 
should resemble Gould's prediction of the DF if diffusion 
is the dominant gravitational process in the solar system, 
which would yield a bound WIMP number density not 
significantly larger than predicted in this work (although 
the fact that many of these WIMPs have low speeds has 
disproportionate consequences at neutrino telescopes). 

However, diffusion is not necessarily the dominant 
gravitational process, at least in terms of the size of 
the DF. Neither Gould nor Lundbcrg & Edsjo 's treat- 
ment of the gravitational processes in the solar system 
incorporate resonances (Lundberg & Edsjo 's treatment 
of resonances extends only to the timescale for WIMPs 
to enter the Sun; the treatment of WIMP-planet encoun- 
ters is purely local). While Gladman ct al. [37] find for 
their sample of NEOs, resonances do not change the semi- 
major axis distribution relative to diffusion for a < 2 AU, 
they can affect the other orbital parameters. Since the 
DF is sensitive to how WIMPs encounter the Earth, it is 
important to understand the differences between purely 
diffusive orbital evolution and the evolution when reso- 
nances are present. Moreover, if the resonances can shield 
WIMPs from close encounters with planets (for example;, 
the Kozai resonance assures that either the eccentricity 
or inclination of an orbit is high, making encounters with 
planets only occur at high speeds), they can extend the 
lifetimes of certain classes of orbits. It is important to 
determine if resonances yield a significant long lifetime 
tail in the Earth-crossing WIMP distribution. 

In conclusion, in a more realistic solar system, we ex- 
pect that prograde WIMPs will have a higher density at 
the Earth than retrograde WIMPs, and that the den- 
sity of resonance-sticking Jupiter-crossing WIMPs will 
be reduced. We suspect that the DF of gravitationally 
bound WIMPs lies somewhere between the predictions 
of Gould and Lundberg & Edsjo in the absence of res- 
onances, depending on the details of scattering in the 
Sun, unless Galactic gravitational fields affect the WIMP 
DF as strongly as suggested by our crude treatment in 
Section III. Since resonances are important in the or- 
bital evolution of NEOs and in the population of solar 
captured WIMPs described in Paper I, we suspect that 
they will also be important for gravitationally captured 
WIMPs whose orbits become interior to Jupiter's. How- 
ever, orbits will need to be integrated in the full solar 
system in order to evaulate the effects of these resonances 
on the distribution function. 



B. Future Simulations 

We would like to test these hypotheses with simula- 
tions of WIMPs in a more realistic solar system. How- 
ever, our experiences with simulating orbits in a toy solar 
system have highlighted some potential difficulties in go- 
ing to a more complicated solar system. As is the case 
for simulating the solar-captured WIMP population, the 
subject of Paper I, the primary problem will be to sim- 



ulate a statistically significant number of orbits in finite 
computing time. In the toy solar system, we required 
~ 2 X 10^ CPU-hours for the integration of ~ 10^° WIMP 
orbits, of which only ^ 3 x 10^ became bound, and only 
~ 6000 of which contributed to the DF at the Earth. 
There were only 5 resonance-sticking WIMPs with life- 
times longer than 10 Myr. Since orbits in the real solar 
system display a much richer phenomenology than in the 
toy solar system, the number of WIMPs simidated in the 
toy solar system is insufficient to sample the spectrum 
of behavior in a more realistic solar system; in fact, we 
barely simulated enough WIMPs to find the resonance- 
sticking phenomenon in the toy solar system simulations. 

Just as we did in Paper I (Section VIIB), we propose a 
few techniques for exploring solar system phenomenology 
in finite computing time. First, as in Paper I, we recom- 
mend a series of intermediate simulations leading up to 
a full solar system simiilation. Perhaps the zero-order 
simulation in this series would include a more realistic 
treatment of Galactic gravitational fields in another toy 
solar system simulation, and a sufficiently large number 
of particles to determine the resonance-sticking WIMP 
DF. The next simulation, in order to explore phenomena 
associated in systems with multiple planets, would in- 
clude just three planets: Jupiter, Earth, and Venus. The 
planets would initially be put on circular, coplanar orbits 
about the Sun to highlight diffusion phenomena as well 
as mean-motion and Kozai resonances. Since the gravi- 
tational cross section scales as the square of the planet 
mass, we propose scaling up the masses of the Earth and 
Venus. The masses of the planets could be reduced in 
follow-up simulations. A possible next step would be to 
put the limited set of planets on more realistic orbits in 
order to explore changes to diffusion as well as the effects 
of secular resonances. 

Second, we propose using more clever choices for the 
initial conditions of the simulations. While we weighted 
the initial conditions towards low energy WIMPs and re- 
stricted the range of angular momenta, still only a small 
fraction of WIMPs were captured in the solar system in 
order to boost the yield of bound WIMPs. We sam- 
pled WIMPs up to energies of ~ 0.5(44 km s~^)^. 
Since the angular momenta were quite high in the High 
Perihelion run, the only WIMPs captured in the solar 
system had energies several orders of magnitude below 
this upper limit. In the future, we recommend restrict- 
ing High Perihelion-like simulations to lower energies. In 
general, though, since WIMP orbits in the solar system 
are chaotic and the gravitational cross sections are small, 
it is difficult to further improve the yield of bound. Earth- 
crossing orbits by the choice of the weighting and range 
of E and for the incoming halo WIMPs. 

However, if it appears that the long lifetime tail of 
the DF in the toy solar system and additional capture of 
WIMPs from the halo by the outer planets are unimpor- 
tant, initial conditions for the simulations may be drawn 
from the toy solar system DF. This is an attractive choice 
because the cross section for Jupiter to capture a halo 
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WIMP is small (requiring us to simulate ~ 10^" halo 
WIMPs to obtain a sample of ~ 10^ bound WIMPs); 
by sampling the toy solar system DF, a far higher frac- 
tion of the simulated WIMPs should contribute to the 
full solar system WIMP DF at the Earth. The reason 
to think this might work is because the typical diffusion 
timescales for the other planets in the solar system are 
far longer than for Jupiter. Therefore, for a short period 
of time after the birth of the solar system, the toy solar 
system DF would be an accurate representation of the 
full solar system DF. This technique is used in Lundberg 
and Edsjo [31] to speed up their orbital diffusion cal- 
culations in a Jupiter- Venus-Earth solar system. How- 
ever, there are two potential drawbacks to this approach. 
First, if the long lifetime tail in the toy solar system is im- 
portant, the characteristic timescales of processes in the 
solar system related to planets other than Jupiter may 
approach the characteristic timcscalc in the DF, meaning 
that the long lifetime tail could be significantly altered 
by the other planets. Second, if secular resonances act on 
short timescales to pump WIMPs into the Sun, as sug- 
gested in simulations of NEO orbits [35, 37], the toy solar 
system DF will not be an accurate representation of the 
full solar system DF even on very short timescales. Care- 
ful tests need to be performed to evaluate the possibility 
of this approach to the initial conditions. 

In summary, we emphasize the need for a series of sim- 
ulations, culminating in a full solar system simulation, in 
order to understand the importance of different phenom- 
ena in the solar system. The choice of initial conditions 
is important to maximize statistics, but we are limited 
by the chaotic nature of the system. 



C. The Halo DF 

Gravitational capture is most effective for halo WIMPs 
with speeds ^ 1 km s'~^ outside the potential well of the 
Sun. Therefore, it is important to determine the uncer- 
tainty on the halo DF at such small heliocentric speeds. 
We briefly detail three ways in which deviations from our 
fiducial halo WIMP model may affect the gravitational 
capture rate N of WIMPs in the halo, and hence, the 
bound WIMP DF. Since we discussed anisotropic halo 
speed distributions in Section III, we omit further refer- 
ence to anisotropy here. 

1. Uncertainty in Vq, a, and p-^: The halo model (Eq. 
25) used to set up the initial conditions (Section II B) 
and to derive the WIMP DF from simulation outputs 
(Appendix A) is a single-variate Gaussian, non-rotating 
in a Galactocentric frame, and assumes that the speed 
of the Local Standard of Rest (LSR) vq = 220 km 
and the one-dimensional halo WIMP velocity dispersion 
is (7 = vq/\^. The uncertainty in the speed of the LSR 
is ~ 10% [60]. While the fiducial value of Vq is what 
was used in our fiducial halo model, recent astrometry of 
masers in star-forming regions in the spiral arms of the 
Milky Way suggests Vq « 250 km s~^ [61]. To determine 



what this uncertainty in vq implies for the uncertainty 
in the capture rate of WIMPs to the solar system, we 
find that the angle-averaged halo DF outside the sphere 
of influence of the Sun (average over Eq. 26) 

for the heliocentric speed <C w©. If cr is fixed, a 10% 
uncertainty in vq implies a ~ 20% uncertainty in the 
low speed halo WIMP DF, and hence, a 20% uncertainty 
in the capture rate of WIMPs from the halo. We note 
that we have neglected solar motion in our fiducial model, 
which is also of order 10% of the speed of the Local Stan- 
dard of Rest [62]. 

However, we set a = Vq/ \pi in the halo DF. Therefore, 
if we include the uncertainty of a through the uncertainty 
in U0, we find that a 10% uncertainty in Vq yields a 30% 
uncertainty in the DF. 

There is further uncertainty in a. We used a = 
since this is the one-dimensional velocity dispersion for an 
isothermal halo DF for a power-law density distribution 
n^{r) oc r~^. Since halos are approximately described by 
an NFW profile, 

n^(r) (X (r/r,)-i(l+r/r,)-2, (51) 

where is a scale radius, n^(r) oc corresponds to 
the transition between the inner cusp {n^{r) (x r~^) and 
the outer halo {n^{r) cx r^^) [63]. Simulations suggest 
that the solar circle should be near this transition zone, 
although only small changes in the position of the Sun 
with respect to the transition radius can yield different 
relationships between the circular speed and the velocity 
dispersions (see Fig. 11 in [64]). In general, a power- 
law density distribution ny^{r) oc r^^^ yields a velocity 
dispersion a = Vq/ for an isothermal DF, so if the Sun 
is well within the transition, a = Vq^ while u = Vq/\/3 
if outside (see Appendix A of [65]). 

If there is no uncertainty in vq and it is fixed to the 
fiducial value, the uncertainty in the relation between 
a and vq yields an uncertainty in the low speed halo 
WIMP DF of order 40%. Thus, the total uncertainty 
in the gravitational capture rate of WIMPs to the solar 
system due to uncertainties in a and Vq is of order 50% 
assuming that the WIMPs have an isotropic MaxweUian 
velocity distribution. 

Simulations show that a multi-variate Gaussian is a 
good fit the macroscopic velocity distribution near the 
solar system at the speeds relevant for capture in the solar 
system, while we have used a single-variate Gaussian [26, 
64, 68] . However, the velocity dispersion in each direction 
is within ^ 10 — 20% of the one-dimensional velocity 
dispersion used in our simulations. Therefore, we believe 
that the uncertainty in f3 yields greater uncertainty in 
the capture rate than treating the halo as a single-variate 
instead of a multi-variate Gaussian. 

Finally, we consider uncertainties in the local WIMP 
density p^. Dark matter-only N-body simulations show 
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that the density at any radius from the Galactic center 
should be within ~ 10 — 20% of its mean value in the 
shell. However, the observed uncertainty in the local dark 
matter density is much higher because baryons dominate 
the local potential and density, and the distribution of 
baryons in the Galaxy is somewhat uncertain [69, 70]. 
There is approximately a factor of two uncertainty in 
the local dark matter density owing to the uncertainty 
in the baryon distribution [25]. This yields a factor of 
^ 4 uncertainty in the annihilation rate of WlMPs in 
the Earth. 

2. Halo substructure: Finally, we consider the issue of 
dark matter substructure. Substructure will have a small 
effect on the DF of gravitationally captured WIMPs for 
the same reasons why it has a negligible effect on the DF 
of solar captured WIMPs (discussed in Paper 1). Briefly, 
the probability that the solar system is in a dense subhalo 
is ^ 10~^ at any given time [68]. The equilibrium time 
of the gravitationally captured WIMP DF is sufficiently 
large that the solar system should have passed through 
many clumps in the equilibrium time. The velocity dis- 
tribution of subhalos is somewhat biased, at the level of 
the velocity dispersion of subhalos being ~ 50% higher 
than the smooth component [71, 72]. This yields order 
unity or less changes in the capture rate relative to the 
capture rate if the velocity distribution of subhalos were 
unbiased with respect to the smooth component. Thus, 
the capture rate of WIMPs to the solar system averaged 
over the DF equilibrium time (N) ~ N, where N is the 
capture rate of the smooth component in the absence of 
substructure. 

Thus, we find that the presence of substructures in 
the halo and the uncertainty in vq, a, and p-^ of the 
smooth halo WIMP DF only change the overall capture 
rate of WIMPs in the solar system by at most of order 
unity. This will yield a factor of a few uncertainty in the 
annihilation rate of WIMPs in the Earth, although it is 
unlikely that the rate will be high enough to produce an 
observable signal in IceCube. 

3. Macroscopic dark matter structure: In hydrody- 
namic N-body simulations of the evolution of disk galax- 
ies, a second dark matter structure has been found. This 
is a "thick disk" of dark matter, whose properties mimic 
the thick stellar disk, and is a result of satellite galaxies 
being preferentially dragged into the disk plane as they 
merge with disk galaxies [27, 73, 74] . The thick disk has 
of order the same density at the solar circle as the more 
spherical dark matter halo, but is oblate and rotates in 
the same sense as the Sun about the Galactic Center. 
The dark matter in the thick disk typically has a smaller 
speed with respect to the Sun and a smaller velocity dis- 
persion than halo dark matter. Both effects tend to boost 
the capture rate of WIMPs in the solar system. The im- 
plications of this structure on the bound WIMP DF and 
event rates in neutrino telescopes are described in Bruch 
et al. [75]. 



VII. CONCLUSION 

The main conclusions of this work are: 

• We found the DF of WIMPs gravitationally bound 
to a toy model solar system consisting of Jupiter 
on a circular orbit about the Sun. 

• While the DF matches the detailed balance pre- 
diction for geocentric speeds ?; < 30 km s~^, the 
fit to that prediction is poor for larger speeds. 
This discrepancy is largely due to the difficulty 
in populating retrograde. Earth-crossing orbits. 
Wc find a small population of long-lived Earth- 
crossing resonance-sticking orbits, contributing ~ 
15% to the number density of gravitationally bound 
WIMPs at the Earth. 

• The DF of gravitationally bound WIMPs is largely 
insensitive to the details of WIMP-nucleus interac- 
tions in the Sun, confirming the findings of Lund- 
berg and Edsjo [31]. 

• If external Galactic gravitational fields become im- 
portant for distances r J5J 1000 AU from the Sun, 
the DF could be significantly smaller than shown 
in Fig. 7, by up to a factor of three. Future simula- 
tions of WIMPs in the solar system should include 
more realistic treatments of external gravitational 
fields. 

• The DF of gravitationally captured WIMPs is of 
order the same size as the largest DF of solar 
captured WIMPs consistent with limits on the 
WIMP-nucleon cross section (Paper I). Since the 
DF of gravitationally captured WIMPs is larger at 
lower speeds, it dominates the annihilation rate of 
WIMPs in the Earth for > 100 GeV. Since the 
DF of solar captured WIMPs is larger at higher 
geocentric speeds, solar captured WIMPs dominate 
the direct detection signal of bound WIMPs. 

• The direct detection event rate of gravitationally 
captured WIMPs is never more than ^ 0.1% of the 
event rate of halo WIMPs. 

• We find that the annihilation rate in the Earth of 
WIMPs gravitationally captured to the solar sys- 
tem agrees with Gould [33] 's prediction for a toy 
solar system. However, even though gravitationally 
captured WIMPs enhance the event rate of neutri- 
nos from WIMP annihilation in the Earth's core, 
the signal falls far short of threshold for the Ice- 
Cube experiment assuming a standard halo model. 

Acknowledgments 

We thank Scott Tremaine for advising this project, and 
Aldo Serenelli and Carlos Penya-Garay for providing ta- 
bles of isotope abundances in the Sun. We acknowledge 



22 



financial support from NASA grants NNG04GL47G and 
NNX08AH24G and from the Gordon and Betty Moore 
Foundation. The simulations were performed using com- 
puting resources at Princeton University supported by 
the Department of Astrophysical Sciences (NSF AST- 
0216105), the Department of Physics, and the TIGRESS 
High Performance Computing Center. 



APPENDIX A: ESTIMATING DISTRIBUTION 
FUNCTIONS 

In this Appendix, we describe how to estimate the 
WIMP DF at the Earth from the simulations, taking into 

account the weighting of initial conditions relative to the 
flux of halo WIMPs through the solar system. Instead of 
estimating the WIMP DF at each point along the Earth's 
orbit, we find the DF averaged along the Earth's orbit, 
and within a height Zc ^ a® of the reference plane (as- 
suming the orbits of the Earth and Jupiter are coplanar). 
By averaging over a small region along the Earth's orbit, 
we improve the DF statistics; by choosing Zc to be small, 
the averaged DF should not be contaminated by gradi- 
ents in the DF as a function of height above the reference 
plane. 

To construct the averaged DF, we record the phase 
space coordinates each time a WIMP passes through a 
cylindrical shell of radius a® and height Izc centered on 
the reference plane. Thus, in effect, we record the unnor- 
malized WIMP flux as a function of time. To find the 
DF from these data, we must weight the flux with respect 
to the initial conditions of the simulation, and relate the 
flux to the DF. 

The flux F{y) is related to DF as follows. The number 
of WIMPs passing outward through a wall of size 
oriented in direction in time 5t is 



6N = — dvSASt, 
dv 



(Al) 



where v is the WIMP speed in a frame in which the wall 
is fixed. Equivalently, the WIMPs that pass through the 
wall in time St inhabit a prism of base 6 A, long side vSt, 
and height Stv-6A with a number density /(v)dv. There- 
fore, in terms of the DF /(v), the number of WIMPs 
passing outward through the wall is 



SN = f{v)dv6A6tiv-dA) 
= f{v)v±dv6A6t, 



(A2) 
(A3) 



where v± is the component of WIMP speed perpendicular 
to the wall with outward normal 6 A. Remembering that 
the DF is always positive, we equate Eqs. (Al) and (A3), 
we find 



/(v) 



di^(v)/dv 



(A4) 



where v± is component of the WIMP velocity perpendic- 
ular to the surface of the cylinder. 



Therefore, if we find the differential fiux dF/dv, we can 
quickly determine the WIMP DF. First, we must weight 
the WIMPs in the flux according to their initial condi- 
tions. In Section II B, we weighted the WIMP initial 
conditions relative to the flux of halo WIMPs through 
the solar system by a factor of 



W{E) = 



j,max ^ 



(A5) 



assuming an isotropic halo WIMP velocity distribution. 
Therefore, we weight the WIMPs in the flux by 



w{E) = W-'^{E) 



J^iEr, 



-(A6) 



This would be the proper weighting if the typical time for 
the WIMP DF to reach equilibrium were longer than the 
time for the Sun to circle the Galactic center, ^ 200 Myr, 
since the time averaged DF of Eq. (26) is nearly isotropic. 
However, in Section III, the equilibrium timescale is less 
than this. Thus, the particles should be weighted by 
their initial position and velocity relative to the orbital 
planet of the solar system and not just their initial speed 
(or energy). For now, though, we treat the halo DF as 
isotropic. 

The estimated average WIMP flux through the cylin- 
der, properly normalized with respect to the initial par- 
ticle distribution, is given by 



dF(v, t) N 



dv 2Tra(nZ, 



^ ^u;(i;«)(5(3)(v-v«^)e(t-i„^) 



a=l/3=l 



/Y.w{E^). (A7) 



a=l 



Here, N is the rate at which halo WIMPs in the energy 
and perihelion range considered in a simulation enter the 
solar system, Np is the total number of WIMP simulated 
in each run, a labels a particular WIMP, and [3 denotes a 
particular passage of WIMP a through the cylinder. Vq,^ 
is the velocity of the WIMP a during passage (3, and tdp 
is the time of that passage since the birth of the solar 
system. The denominator in Eq. (A7) normalizes the 
flux. To find the estimated bound WIMP DF, we insert 
Eq. (A7) into Eq. (A4). 

To find the DF, we set z^ = lO-^a^ w 23.5i?e. The 
DF does not change within confidence limits for smaller 
values of Zc- To estimate uncertainties the DF, we employ 
bootstrap resampling, drawing Np WIMPs with replace- 
ment from the initial WIMP sample. For each simulation 
run, we do this 500 times. To find the total DF of WIMPs 
bound to the solar system by gravitational capture, we 
sum the DFs from the Regular and High Perihelion runs 
and add the errors in quadrature, since the simulation 
runs are independent of each other. 
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